home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Skunkware 5
/
Skunkware 5.iso
/
lib
/
calc
/
chrem.cal
< prev
next >
Wrap
Text File
|
1995-07-17
|
4KB
|
182 lines
/*
* chrem - Chinese remainder theorem/problem solver
*
* When possible, chrem finds solutions for x of a set of congruences
* of the form:
*
* x = r1 (mod m1)
* x = r2 (mod m2)
* ...
*
* where the residues r1, r2, ... and the moduli m1, m2, ... are
* given integers. The Chinese remainder theorem states that if
* m1, m2, ... are relatively prime in pairs, the above congruences
* have a unique solution modulo m1 * m2 * ... If m1, m2, ...
* are not relatively prime in pairs, it is possible that no solution
* exists. If solutions exist, the general solution is expressible as:
*
* x = r (mod m)
*
* where m = lcm(m1,m2,...), and if m > 0, 0 <= r < m. This
* solution may be interpreted as:
*
* x = r + k * m [[NOTE 1]]
*
* where k is an arbitrary integer.
*
***
*
* usage:
*
* chrem(r1,m1 [,r2,m2, ...])
*
* r1, r2, ... remainder integers or null values
* m1, m2, ... moduli integers
*
* chrem(r_list, [m_list])
*
* r_list list (r1,r2, ...)
* m_list list (m1,m2, ...)
*
* If m_list is omitted, then 'defaultmlist' is used.
* This default list is a global value that may be changed
* by the user. Initially it is the first 8 primes.
*
* If a remainder is null(), then the corresponding congruence is
* ignored. This is useful when working with a fixed list of moduli.
*
* If there are more remainders than moduli, then the later moduli are
* ignored.
*
* The moduli may be any integers, not necessarily relatively prime in
* pairs (as required for the Chinese remainder theorem). Any moduli
* may be zero; x = r (mod 0) has the meaning of x = r.
*
* returns:
*
* If args were integer pairs:
*
* r ('r' is defined above, see [[NOTE 1]])
*
* If 1 or 2 list args were given:
*
* (r, m) ('r' and 'm' are defined above, see [[NOTE 1]])
*
* NOTE: In all cases, null() is returned if there is no solution.
*
***
*
* This function may be used to solve the following historical problems:
*
* Sun-Tsu, 1st century A.D.
*
* To find a number for which the reminders after division by 3, 5, 7
* are 2, 3, 2, respectively:
*
* chrem(2,3,3,5,2,7) ---> 23
*
* Fibonacci, 13th century A.D.
*
* To find a number divisible by 7 which leaves remainder 1 when
* divided by 2, 3, 4, 5, or 6:
*
*
* chrem(list(0,1,1,1,1,1),list(7,2,3,4,5,6)) ---> (301,420)
*
* i.e., any value that is 301 mod 420.
*
* Written by: Ernest W Bowen <ernie@neumann.une.edu.au>
* Interface by: Landon Curt Noll <chongo@toad.com>
*/
static defaultmlist = list(2,3,5,7,11,13,17,19); /* The first eight primes */
define chrem()
{
local argc; /* number of args given */
local rlist; /* reminder list - ri */
local mlist; /* modulus list - mi */
local list_args; /* true => args given are lists, not r1,m1, ... */
local m,z,r,y,d,t,x,u,i;
/*
* parse args
*/
argc = param(0);
if (argc == 0) {
quit "usage: chrem(r1,m1 [,r2,m2 ...]) or chrem(r_list, m_list)";
}
list_args = islist(param(1));
if (list_args) {
rlist = param(1);
mlist = (argc == 1) ? defaultmlist : param(2);
if (size(rlist) > size(mlist)) {
quit "too many residues";
}
} else {
if (argc % 2 == 1) {
quit "odd number integers given";
}
rlist = list();
mlist = list();
for (i=1; i <= argc; i+=2) {
push(rlist, param(i));
push(mlist, param(i+1));
}
}
/*
* solve the problem found in rlist & mlist
*/
m = 1;
z = 0;
while (size(rlist)) {
r=pop(rlist);
y=abs(pop(mlist));
if (r==null())
continue;
if (m) {
if (y) {
d = t = z - r;
m = lcm(x=m, y);
while (d % y) {
u = x;
x %= y;
swap(x,y);
if (y==0)
return;
z += (t *= -u/y);
}
} else {
if ((r % m) != (z % m))
return;
else {
m = 0;
z = r;
}
}
} else if (((y) && (r % y != z % y)) || (r != z))
return;
}
if (m) {
z %= m;
if (z < 0)
z += m;
}
/*
* return information as required
*/
if (list_args) {
return list(z,m);
} else {
return z;
}
}
global lib_debug;
if (lib_debug >= 0) {
print "chrem(r1,m1 [,r2,m2 ...]) defined";
print "chrem(rlist [,mlist]) defined";
}